% Generate the confidence intervals reported in Table 3.
% Execute separately for lambda = 0.3 and lambda = 0.

clear
clc

% Designate working directories
dir = '.../Replication package/';
cd(dir);
addpath(genpath(strcat(dir,"Code/src")))
addpath(genpath(strcat(dir,"Datasets")))
filename    = strcat(dir,'Table_3.txt');
fileID      = fopen(filename, 'w');

% Set seed number
seed = 4534546 ;

bnd          = 5 ;
lambda       = 0.3;
BR_P         = 1000 ;
BR_G         = 5000 ;
sz_Gamma     = 21^4 + 1 ;
theta_grid   = -0.5:0.01:0.5 ;
splits       = 1000 ;
swarmsize    = 200 ;
tol          = 1e-2 ;
alpha_level  = 0.05 ;  % 95% CI is computed

% Load the dataset
data = table2array(readtable('USAY.txt'));
data = data(2:end, 1:8);
T    = size(data, 1) - 2;
Y    = data(3:end, 6);
X    = data(3:end, 8);
Z    = [data(1:(end-2), 3:6)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use PSO algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
eps_constant = 1e-2;
pso_nmax     = 50;

rng(seed,'twister');

ts    = zeros(length(theta_grid),1);
ts_bs = zeros(BR_P, length(theta_grid),1);
zeta  = normrnd(0, 1, [T, BR_P]);

for i = 1:length(theta_grid)
    datau = Y-mean(Y) - theta_grid(i)*(X-mean(X));
    ts(i, :) = pen_Bierens_test(datau,Z,bnd,lambda,eps_constant,pso_nmax,'swarmsize',swarmsize,'dm',1);
    for b = 1:BR_P
        ts_bs(b, i, :) = pen_Bierens_test(zeta(:, b).*datau,Z,bnd,lambda,eps_constant,pso_nmax,'swarmsize',swarmsize,'dm',1);
    end
end

cv          = reshape(quantile(ts_bs, 1-alpha_level, 1), [length(theta_grid) 1]);
fail_to_rej = (ts <= cv);

CI_PSO = theta_grid(fail_to_rej);

fprintf(fileID, '%.2f%% confidence interval computed using PSO : [%.2f, %.2f]\n', 100*(1-alpha_level), min(CI_PSO), max(CI_PSO));
fprintf(fileID, '---------------------------------------------------------------------------\n');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use grid search method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rng(seed,'twister');

out_GS = pen_Bierens_test_GS(Y-mean(Y)-(X-mean(X))*theta_grid,Z,bnd,lambda,'bR',BR_G,'sz_Gamma',sz_Gamma,'splits',splits,'seed',seed,'alpha_level',alpha_level,'dm',1);

CI_GS = theta_grid(out_GS.rejection == 0);

fprintf(fileID, '%.2f%% confidence interval computed using GS : [%.2f, %.2f]\n', 100*(1-alpha_level), min(CI_GS), max(CI_GS));
fclose(fileID);